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Abstract 

We study the Ising model on the triangular lattice with nearest-neighbor couplings K nn , next- 
nearest-neighbor couplings K nnn > 0, and a magnetic field H. This work is done by means of 
finite-size scaling of numerical results of transfer matrix calculations, and Monte Carlo simulations. 
We determine the phase diagram and confirm the character of the critical manifolds. The emphasis 
of this work is on the antiferromagnetic case K nn < 0, but we also explore the ferromagnetic regime 
Knn > for H = 0. For K nn < and H = we locate a critical phase presumably covering the 
whole range — oo < K nn < 0. For K nn < 0, H ^ we locate a plane of phase transitions containing 
a line of tricritical three-state Potts transitions. In the limit H — * oo this line leads to a tricritical 
model of hard hexagons with an attractive next-nearest-neighbor potential. 

PACS numbers: 05.50.+q, 64.60.Ak, 64.60. Cn, 64.60.Fr 
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I. INTRODUCTION 



The Ising model on the triangular lattice with nearest-neighbor couplings K nn , next- 
nearest-neighbor couplings K nnn , and a magnetic field H, is defined by the reduced Hamil- 
tonian 

H/k B T = -K nn ^ S i S j ~ K mm ^ SkSl - H ^2 S m, (1) 
(nn) [nnn] m 

where = ±1, and (nn) and [nnn] indicate summations over all pairs of nearest neighbors 
and of next-nearest neighbors, respectively, as illustrated in Fig. ^ 

This model, in particular the antiferromagnetic model (K nn < 0), displays interesting 
behavior. For K nnn = 0, H = the model has been solved exactly A ferromagnetic 
transition occurs at K nn = ln(3)/4. An antiferromagnetic (K nn < 0) mirror image of this 
transition is absent. This is related to the fact that the triangular lattice is not bipartite. 
However, at zero temperature, i.e., for K^a — > — oo, the model displays a critical phase 
with algebraically decaying correlations |2j. This zero-temperature model can be exactly 
mapped on a solid-on-solid (SOS) model |3[. Under renormalization, it is assumed to map 
on the Gaussian model 4] and on the related Coulomb gas [5]. The coupling constant g R 
of the Coulomb gas can thus be obtained exactly as gR = 2 so that a number of critical 
exponents can be calculated. The Ising temperature T oc —K m l appears to be relevant: the 
critical state is destroyed for all T > 0. Commensurate-incommensurate transitions occur 
when finite differences between the infinite nearest-neighbor couplings in the three lattice 
directions are introduced HQ]. 

Next we consider the case of if = and K nnn ^ 0. The mapping on the SOS model 
(and we may also assume this for the Coulomb gas) is still valid for K nn — ► — oo but, in the 
absence of an exact solution, g R is no longer exactly known. It has, however, been deduced 



4] 



that g R is an increasing function of K nnn . The Coulomb gas analysis predicts that, for 
sufficiently large gR, the Ising temperature becomes irrelevant, so that the algebraic phase 
extends to nonzero temperatures. This analysis also predicts that for even larger qr a phase 
transition to a flat SOS phase occurs, both at zero and at nonzero temperatures. 

Somewhat earlier, part of this scenario had already been described by Landau [6]. Via 
the lattice-gas representation of Eq. (JTj), he used the connection with the XY model in the 
presence of a six-state clock-like perturbation, made earlier by Domany et al. 0. He could 
thus make use of their results Q for this model which allow for the existence of a critical, 



XY-like phase in a nonzero range K nn > 0. Furthermore, Landau |6( used the Monte Carlo 
method to verify the existence and nonuniversal character of this critical phase for the case 
of a fixed ratio K nnn /K nn = — 1. 

Another tool to study the model with nonzero next-nearest-neighbor couplings K nnn is 
provided by the transfer-matrix technique. A simplification has been used in the latter 
approach: K nnn was taken to be nonzero only for four out of the six next-nearest neighbors 
P. Q|. This leads to a substantial simplification of the transfer matrix calculations, 
but the resulting system lacks isotropy, so that applications of conformal mappings become 
difficult. On this basis, limited evidence ljj for the existence of the critical phase was 
reported; the limitation of this evidence is obviously related to the lack of sixfold symmetry. 

Next we consider the consequences of a nonzero field H > 0. On the basis of the relation 
with the Coulomb gas it has been derived jj| that, for K nn — > — oo and K nnn = 0, the 
magnetic field H is irrelevant: the critical state is not destroyed by a sufficiently small field 
H ^ 0. However, the magnetic field tends to increase the Coulomb gas coupling constant 
gn. The field will become marginally relevant at Qr = 9/4 and a transition of the Kosterlitz- 
Thouless (KT) type or, in this context more appropriate, of the roughening type is thus 
expected. This transition separates the critical phase from a long-range ordered phase, 
where the majority of the minus-spins have condensed on one of the three sublattices of 



12j by means of numerical 



the triangular lattice. This prediction has been confirmed 111 
methods. The long-range ordered phase extends to nonzero temperature T > and is 
separated from the disordered phase by a line of phase transitions in the (H, T) plane that 
belongs to the three-state Potts universality class [3, Q, 13 3, Q • 

Since the Ising model in a field can be mapped on a vertex model, and the critical 
manifolds of solvable vertex models are described by the zeroes of simple polynomials in the 
vertex weights it may be assumed that also for the triangular lattice the critical line 
in the (H, T) is described by such a polynomial. This assumption was recently refuted by 
Qian et al. [l2|. The shape of the critical line, as deduced from this assumption, was found 
to be inconsistent with the numerical evidence. They also found that the renormalization 
ideas originally outlined by Nienhuis et al. [4] could be applied to predict the shape of the 
critical line in the (H, T) plane for small T. This shape was found to be consistent with 
their numerical data for the critical line. 

The aforementioned three-state Potts-like critical line is naturally part of a critical surface 
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extending to nonzero K nnn . The more involved problem to find the phase diagram in the 
three-parameter (if, K nnn , K nn ) space has already been partly explored. On the basis of 
renormalization arguments, Nienhuis et al. j^J obtained information about the shape of the 
critical surface in the limit H — > 0. Landau p performed Monte Carlo simulations for a 
fixed ratio K nnn /K nn = — 1. He determined the line of phase transitions as a function of H 
and noted that the three-state Potts character along this line changes at a tricritical point 
beyond which the transition turns first-order. 

In this work we verify the predictions in Ref. U and determine the critical values of K nnn 
corresponding to several relevant values of the Coulomb gas coupling constant gn, both for 
finite and infinite K nn . We verify the character of the predicted critical phase at H = 0. We 
also study the critical phenomena associated with the introduction of a nonzero magnetic 
field and explore the full three-parameter phase diagram for K nnn > 0. 

This paper is organized as follows. In Sec. [HJ we summarize our numerical methods which 
include Monte Carlo algorithhms and the construction of a transfer matrix. We define the 
observables that will be the subject of our numerical analysis. The study of the phase 
transitions of the triangular Ising model with nearest- and next-nearest-neighbor couplings 
in a zero field is presented in Sec. IIHI and in Sec. II VI we describe our results for a nonzero 
magnetic field; we conclude with a discussion in Sec. 



II. NUMERICAL METHODS 
A. Transfer-matrix calculations 



Most of the the transfer-matrix calculations were performed for T > so that we had to 
use a binary representation for the Ising spins, leading to a transfer matrix of size 2 L x 2 L for a 
system with finite size L. For T = one can use a simplified transfer matrix of a smaller size 



111] . We define the spin lattice on the surface of a cylinder, and take the transfer direction 



perpendicular to a set of nearest-neighbor edges. The lattice is divided into three sublattices 
denoted as 1, 2 and 3, respectively, as shown in Fig. ^ Nearest-neighbor interactions occur 
only between different sublattices and next-nearest-neighbor interactions occur within the 
same sublattice. 

To enable calculations for system as large as possible, a sparse matrix decomposition 
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has been used. This leads to a very significant reduction of the required computer time 



and memory. The transfer matrices are defined in |11|, [12j for the nearest-neighbor model. 
Here we modify the transfer matrix to include all next-nearest-neighbor interactions. This 
makes it necessary to code two (instead of one) layers of spins as the transfer matrix index. 
Finite-size calculations with L multiples of 6 up to L = 24 were performed. The maximum 
finite size L = 24 corresponds to a cylinder with a circumference of only 12 nearest-neighbor 
bonds. 

The magnetic correlation function along the coordinate r in the length direction of the 
cylinder is defined as: 

9m(r) = (s s r ). (2) 

At large r, this correlation function decays exponentially with a characteristic length scale 
£ m that depends on K nn , K nnn , H, and L: 

g m (r) OC e -r/U(Knn,K nnn ,H,L) (3) 

which can be calculated from the largest two eigenvalues Ao and Ai of the transfer matrix: 

C^nn, K nnn , H,L) = ^= ln(VAi), (4) 

where the factor 2y/3 is a geometric factor for two layers of spins. For the calculation of £ m , 
we make use of the symmetry of the eigenvectors associated with Ao and Ai. The leading 
eigenvector (for Ao) is invariant under a spatial inversion. In contrast, the second eigenvector 
is antisymmetric under inversion. 

The theory of conformal invariance relates £ m on the cylinder with the magnetic 
scaling dimension X m (one half of the magnetic correlation function exponent if). This 
exponent may be estimated as 

X m (K na , K nnn , H, L) = L — — . (5) 

Asymptotically for a critical model with large L we have 

X m (K nn , K nnn , H, L) ~ X m , (6) 

where X m = l/(2gn) in the language of Coulomb gas. This equation allows us to estimate 
X m numerically and thus to obtain evidence about the universality class of the model. Or, if 
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the universality class, and thus X m , are considered known, Eq. (JEJ) and be used to determine 
the critical surface, e.g. to solve for K nnn for given values of K nn , H and L. As a consequence 
of corrections to scaling, the solution will not precisely coincide with the critical point. The 
effects of an irrelevant scaling field u and a small deviation t with respect to the critical 
value of K nn , or K nnn , or H are expressed by 

X m {K nn , K nnn , H, L) = X m + auL Vl + btU» + • • • , (7) 

where a and b are unknown constants, yi is irrelevant exponent and yt is temperature 
exponent. For the solution of the equation X m (K nn , K nnn , H, L) = X m we thus have 
auL Vi + btL Vt 0, so that we expect corrections proportional to L Vi ~ yt in the critical point 
estimates. For instance, for three-state Potts universality one has yt = 6/5 and yi = —4/5 
so that the leading finite-size dependence of the estimated critical points is as L~ 2 . This 
knowledge is helpful for the extrapolation to the actual L = oo critical point. 

In addition to £ m , it is possible to determine a second correlation length £ t describing 
the exponential decay of the energy-energy correlation function. It is associated with a 
third eigenvalue A2 of the transfer matrix with an eigenvector that is symmetric under a 
spatial inversion, just as the one with eigenvalue Ao- The pertinent eigenvalue is thus solved 
by means of orthogonalization with respect to the first eigenvector. In analogy with the 
case of the magnetic correlation length we can use the third eigenvalue A2 to estimate the 
temperature-like scaling dimension X t as 

X t (K nn , K nnn , H, L) = L — , (8) 

where £ t = ^7^1n(Ao/A 2 ). At criticality, it behaves for large L as: 

X t (K Qn ,K Qnn ,H,L)~X t . (9) 

Combining Eqs. © and 0, we can solve for two unknowns simultaneously, using the known 
values of the tricritical three-state Potts model, namely X m = 2/21 and Xt — 2/7. In 
this way, we can estimate the tricritical point (K nnn , K nn ) for a given H. The corrections can 
be argued to be proportional to L m ~ yt2 where y t 2 = 4/7 and yi = —10/7, i.e., the corrections 
decay as L~ 2 . 
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B. Monte Carlo simulations 



Since transfer-matrix calculations are, although highly accurate, restricted to small sys- 
tems, we have also written Monte Carlo algorithms for the present model. To obtain good 
statistical accuracies we included not only a Metropolis algorithm, but also a Wolff and a 
geometric cluster algorithm. Which algorithm is used depends on the location in the phase 
diagram. The Wolff algorithm is applicable in only the case of zero magnetic field. The 
geometric algorithm 19( conserves the magnetization and was therefore used in combina- 
tion with the Metropolis algorithm. This combination was found to work faster than the 
Metropolis method, but the gain in efficiency depends on the position in the three-parameter 
space. 

Several quantities were sampled using these algorithms in order to explore the phase 
diagram. First we define the uniform magnetization as m = L~ 2 ^ k Sk which tends to ±1/3 
in the long-range ordered antiferromagnetic or 'flat' phases, and to zero in the disordered 
(paramagnetic) phase. From its moments we define the magnetic Binder ratio as 

Next, we consider the three-state Potts-like order parameter or, in the language of the present 
Ising model, the three sublattice magnetizations. We denote the magnetization density of 
sublattice i (i — 1, 2, or 3) as mj. On the basis of the staggered magnetizations we write 
the variance of the Potts order parameter as 
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2 — m 2 + m 2 , + m\ — m 1 m 2 — m 2 m 3 — m 3 mi (11) 



s 



and the corresponding dimensionless ratio as 

= (12) 

At criticality, the quantities Q m and Q s scale as a constant plus irrelevant corrections, i.e., 
they converge to a constant as L increases. This property can be used for the determination 
of critical points. 

III. NUMERICAL RESULTS FOR ZERO FIELD 

We restrict this work to ferromagnetic next-nearest-neighbor interactions (K nnn > 0). 
First, we consider the Ising model in a zero field (H = 0), and study the phase diagram in 



(K nnn , K nn ) plane. We distinguish the cases K nn > and K nn < 0. 

A. Results for the ferromagnetic transition (K mi > 0) 

For the Ising model we have X m = 1/8 so that at criticality we expect that asymptotically 
for large L 

X m (K nn , K nnn , 0, L) ~ - (13) 

o 

from which one can estimate critical points e.g. by solving for K nnn at a given value of K nn 
or vice versa. In certain cases, critical points can be determined accurately by extrapolating 
to L = oo . For instance, for K nnn = we obtain the critical value of the nearest-neighbor 
coupling K nn = 0.2746528(10), which is consistent with the exact result K nn = ln(3)/4. The 
results are shown in Fig. EJ 

We also checked that, at the decoupling point (K nn = 0) the critical value of the next- 
nearest-neighbor coupling K nrm equals the exact value ln(3)/4 . The three sublattices, which 
are also triangular lattices, become independent at the decoupling point. 

B. Results for the antiferromagnetic region (K nn < 0) 

At finite K nn < and small K nnn > 0, the model is obviously disordered. As described 
in the Introduction, with increasing K nnn the model is expected to undergo: 1) a Kosterlitz- 
Thouless transition to a critical phase at the point where the Coulomb gas coupling reads 
gjt = 4, and the corresponding value of the magnetic dimension is X m = l/(2gn) = 1/8; 
2) a roughening transition to a flat phase, and the corresponding value of the magnetic 
dimension is thus X m = 1/18 at gu = 9. We have solved K nnn from Eq. (JHJ) for these two 
values of X m , at several fixed values of K nn . The results were extrapolated to L = oo by 
means of three-point fits involving a constant (the estimated value of K nnn ) plus a finite-size 
correction involving a free exponent. The final estimates are included in the phase diagram 
Fig. El They suggest that the two boundaries of the critical phase merge at the decoupling 
point K nn = 0. Our numerical results include a few special points at zero temperature 
(K im — ► — oo). In the renormalization scenario, their meaning is as follows: 

1. For gR = 9/4 we obtain K nnn = 0.0185 (4). This is where the line of roughening 
transitions in the (K nnn , H) plane meets the K nnn axis. 
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2. For gn = 3 we obtain K nnn = 0.0667 (2). This is where the line of three-state Potts 
transitions in the plane perpendicular to the K nnn axis comes in as a straight line with 
a nonzero, finite slope as argued in Ref. [3 

3. For gji = 4 we obtain K nnn = 0.1179 (2). This is where the KT-like line in the 
(K nnn , K nn ) plane meets the K nnn axis. 

4. For gR = 9 we obtain K nnn = 0.226 (2). This is where the line of roughening transitions 
in the (K nnn , K nn ) plane meets the K nnn axis. This point corresponds with an actual 
phase transition on the K nnn axis. We note that, in cases 1 and 3, the K nnn axis 
meets with other lines of phase transitions. However, phase transitions do not occur 
at points 1 and 3 because the critical amplitudes vanish on the K nnn axis. 



C. Shape of the critical lines for small |if nn | 

On the basis of an argument due to van Leeuwen |20j], the scaling behavior of K nn near 
the decoupling point (K nrm = ln(3)/4, K im = 0), is governed by a new critical exponent 
y a = 7/4. This exponent thus determines the shape of the critical lines for small |-K" nn | 
according to 

/In 3 \ 7/4 
K*m oc — K nnn . (14) 
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One can find the critical exponent y a exactly from the known properties of the magnetic 
correlation function of the critical Ising model. The spin-spin correlation behaves as 

g m {r) cc r ~ 2Xm , 

where X m = 1/8 for the 2D Ising model. This also applies to the decoupling point where the 
model decomposes in three independent sublattices. This determines the scaling behavior 
of a four-spin correlation function involving spins in different sublattices in the limit of 

^nnn ~> 

9a(r) = (sooSoiS r0 s r i) = [g m {r)f oc r" 4Xm , (15) 

where soo and soi are nearest-neighbor spins belonging to different sublattices, say sublattices 
1 and 2. The same applies to the pair (s r o, s r i) a t a distance r. Eq. (|15|) describes the energy- 
energy correlation associated with K nn . Its power-law decay is thus expressed by 

g a (r) oc r- 2X % (16) 



where X a is the scaling dimension of the nearest-neighbor energy density. Comparing the 
two Eq. (|T3|) and Eq. (jlBj). we conclude that X a = 2X m = 1/4 and y a = 7/4. 

We verify Eq. (jT3j) by plotting K nn versus [ln(3)/4 — fr nnn ] 7 / 4 for the ferromagnetic critical 
line in Fig. El and for the two lines containing the algebraic phase in the antiferromagnetic 
region in Fig. In all these cases we find approximate linear behavior near the decoupling 
point which confirms the predicted value of y a . 

D. The algebraic phase 

The renormalization scenario predicts that, in the algebraic phase the estimates of X m , 
as obtained from Eq. (|5|). will converge to a if nnn -dependent limit when the finite size L 
increases. However, in the disordered and flat phases, the system will renormalize away 
from the nonuniversal fixed line, and the data for X m are therefore predicted to fan out for 
different values of L. We calculated X m by solving Eq. (J0J) in a suitable range of K nnn at fixed 
values of K nn , namely K nn = — oo, —0.6, —0.4, —0.2 and —0.1. These results confirm the 
renormalization predictions, as illustrated in Figs. El and |H1 Fig.EJshows that, for K nn = — oo 
and H = 0, the data of X m converge to a -ff nnn -dependent constant in a range of K nnn from 
zero to .K'nnn — 0.226 (2) as determined above. This confirms that for H = 0, K nn = — oo the 
system indeed remains critical until K nnn induces a transition to a flat phase. In contrast, 
Fig. IH1 indicates that for nonzero temperature the critical phase starts at a positive value 
of K nnn . Fig. [7| shows the inverse of X m and provides a clearer picture of the transition at 
the large K nnn side. We have numerically calculated the average slopes Sl of the finite-size 
curves in intervals specified in Table I, and fit them as follows: 

S L = S ao + aV" + • • • , (17) 

where is constant, and y c denotes the exponent of the leading finite-size correction. 
Results listed in Table |U indicate that the finite-size dependence of the slopes is governed by 
a negative exponent y c of L, which indicates that the slope Sl converges to a constant for 
L — > oo, as expected in the critical range. 

In order to provide independent confirmation of the algebraic phase, we also used the 
Monte Carlo method. Simulations were done for L x L systems of size L = 24, 36, 48, and 
60. Examples of the results for Q s and Q m are given in Fig. |H] and Fig. El respectively, as a 
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function of K nnn , for K nn = —0.2. These data behave similarly as those for X m , and show 
good apparent convergence to a nonuniversal, i^ nnn -dependent constant in the pertinent 
range. Note that the curves for Q s display intersections near K nnn pa 0.207, and those for 
Q m near K mm pa 0.245, apparently at different sides of the algebraic phase as shown in 
Fig. |21 We interpret these intersections, i.e., solutions of Eq. (JBJ) coinciding for different L, 
as the cancellation of the leading two L-dependent terms. Such terms are likely associated 
with 1) the corrections as naturally associated with irrelevant fields in the algebraic phase; 
and 2) the 'fanning-out' phenomenon mentioned above. It appears that the first types of 
corrections in Q s and Q m are of a different sign. 



IV. RESULTS FOR NONZERO FIELD 

In view of the Ising character of we restrict ourselves to if > without loss of gen- 
erality. The phase diagram without next-nearest-neighbor interactions, i.e., in the (H,K nn ) 
plane has already been determined by Qian et al. |l2j, with special emphasis on the limit 
K Qn — > — oo. In that limit, a roughening-type transition is located Il2| near H = 0.266. 
As mentioned above, the algebraic phase becomes less stable against perturbation by H 
when K nnn increases, and the algebraic phase in the (K nnn , H) plane shrinks to zero at 
Qr = 9/4 which corresponds, as mentioned above, to K nnn = 0.0185. 

The line connecting the two points (K nim , H) =(0,0.266) and (0.0185,0) is a line of rough- 
ening transitions separating the algebraic and the ordered phases. The renormalization de- 
scription implies that this line is a straight line when expressed in the scaling fields. In 
view of the proximity of both numerically determined points, we expect an almost straight 
line in the (K nriri , H) plane. The connection of the three-state Potts transition line and the 
roughening transition point in (H, K nn ) plane has been analytically investigated by Qian et 
al. using renormalization arguments. Their analysis indicates that the roughening transition 
at H = 0.266 is the end point of the Potts transition line in (H, K nn ) plane for T J. 0. Their 
result applies similarly to other points on the line of roughening transitions. We thus believe 
that this whole line serves as a frontier of the Potts critical surface, as well as the part of 
the -fr nn n axis with gn between 9/4 and 4 as determined in Sec. IIII Al 

Since three-state Potts universality implies X m = 2/15 at criticality, we expect that 
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asymptotically for large L 




2 



(18) 
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from which one can estimate critical points by solving for one of the three variables 
(K nn , K nnn , H) for specified values of the other two, and subsequent extrapolation to L = oo. 
We thus calculated critical points on several lines at fixed values of H . The results are shown 
as lines connecting these points in Fig. El I n order to zoom in on the connection of the 
three-state Potts transition surface and the transition lines in the (K nnn , K nn ) plane, we 
have also estimated critical values of H at fixed values of K nn , for a suitably chosen range 
of K nnn . Results for K nn = —0.8, —0.1, —0.15 are included in Fig. E3 They fit well with 
the qualitative predictions for the shape of the critical surface [4] for small H . Furthermore, 
our data for the critical points at K nnn = 0.0667, corresponding with g r = 3, agree with the 
linear behavior as mentioned in Sec. IIII Al 

Our results confirm that, when the next-nearest-neighbor coupling K nnn becomes suffi- 
ciently strong, the transition from the disordered phase to the ordered phase changes char- 
acter at a tricritical line, beyond which the transition turns first-order. We have located the 
tricritical line using transfer-matrix calculations. By solving Eqs. (JBJ) and (JH1) simultaneously 
for K nn and K xmn at specified values of H, we obtain results shown in Table ITT| and included 
in Fig. I n comparison with transfer-matrix calculations involving only X m , the memory 
requirements are somewhat larger. As a consequence only 3 values of L up to 18 could be 
used. But we found that finite-size corrections are relatively small, and we are confident 
that the tricritical line is well determined. 

For sufficiently large fields H, triangles may contain at most one minus-spin and the 
tricritical line approaches a tricritical lattice-gas limit. In this limit the nearest-neighbor 
coupling and the field satisfy a linear relation 



As illustrated in Fig. ^2 the numerical data fit this expression well, except at small H. 
In order to obtain a satisfactory fit to the numerical data for H > 1, we added terms 
proportional to e~ 2H/3 and e ~ 4H/3 to Eq. (|T5jl . This fit yielded C = -0.01481 (5). A similar 
fit without a term proportional to H yielded K nnn = 0.23514 (7) for the tricritical lattice 
gas limit. 



nn 




(19) 
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We have used Monte Carlo simulations to determine the location of the sheet of first- 
order transitions at K nnn = 0.3. We found that, depending on K nn and H, a randomly 
initialized system evolved to a phase either largely magnetized, or resembling one of the 
three ordered Potts states. The threshold values between these two regimes are shown by 
the heavy dashed lines in Fig. They fit smoothly with the results obtained in the critical 
range and for the tricritical line. 



V. DISCUSSION 

We have determined the phase diagram of the model Eq. (JTJ) for K nnn > 0. We locate 
a surface of phase transitions. This surface divides into a three-state Potts-like critical 
sheet and a first-order part. The two parts are separated by a tricritical line. While the 
determination of tricritical line becomes less accurate for small |-K" n n|, our data suggest 
that it spans the whole range — oo < K nn < 0. This is in agreement with the minimal 
renormalization scenario in which the tricritical line is a flow line leading directly from the 
decoupling point to the tricritical fixed point. 

For H — > oo, minus-spins are excluded on nearest-neighbor sites and the the substitution 
o"j = (1 — Sj)/2 reduces the model to a hard-hexagon lattice gas described by the reduced 
Hamiltonian 

5> m , (20) 

(nn) [nnn] m 

where the site variables assume values G{ = 0, 1 and V nn — > oo so that nearest-neighbor 
exclusion applies. The chemical potential of the lattice-gas particles depends on the Ising 
parameters as n — —12K nn — 12K nnn — 2H, and the next-nearest-neighbor potential as 
Vnnn = —4K nnn . For V nnn = this model reduces to Baxter's hard-hexagon lattice gas j^lj. 
According to the analysis presented in Sec. IIVI the tricritical line persists in the lattice-gas 
limit. The Ising parameters C and K nnn determine the tricritical parameters of the lattice- 
2.644 (1) and V„ nn = —0.9406 (3). Our findings may be compared with those of 
Verberkmoes and Nienhuis [22j for a model with V nnn = but including additional smaller 
hexagons. They also report a tricritical point, attributed to an effective attraction between 
the hard hexagons, induced by entropic effects associated with the small hexagons. 

An Ising-like tricritical point is known to occur also in the analogous case of the hard- 



square lattice gas 



23, 



24j . Our result thus confirms that tricriticality is a generic property 
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of hard-core lattice gases with attractive next-nearest-neighbor interactions. 

Since we do not doubt the universality class of the tricritical line, we have not explicitly 
determined its critical exponents. However, we remark that the fast apparent convergence 
of the estimated tricritical points confirms that the values of the Potts tricritical exponents 
X m and X t , as used to solve Eqs. (JBJ) and Q, do indeed apply. 

Renormalization analysis predicts that the uniform magnetic field H is relevant, except 
for a small range 2 < < 9/4. Thus the plane H = qualifies as a possible locus of new 
universality classes, in line with the existence of a critical phase such as predicted by the 
renormalization scenario and confirmed numerically. We finally note that the renormaliza- 
tion equations for the KT transitions imply that the line of KT transitions, as shown in 
Fig. |21 on the left hand boundary of the critical phase, should come in as a straight line 
on the horizontal axis, in contrast with the numerical results which there display a small 
part with a sudden curvature. We believe that this is a finite-size effect, explained by the 
same renormalization equations, which involve the marginally irrelevant temperature field 
parametrizing the line of KT transitions. This scaling field generates slowly converging 
finite-size corrections. This field and its associated finite-size effects vanish at K nn = — oo. 
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TABLE I: Fitted results for the extrapolated average slope Soo ~ dX m /dK miIi in the algebraic 
phase. The last column shows the the exponent y c of finite-size correction. The increase of | S^l with 
-f^nnn corresponds with the narrowing of the algebraic phase when the decoupling point K nn = 
is approached. The intervals of K nnn in which the average slopes are calculated are listed in the 
second column. 



K 


K 


<S*oo 


Vc 


-OO 


0.18 - 0.20 


-0.59 (3) 


-1.1 (2) 


— 0.6 


0.18 - 0.20 


-0.78 (2) 


-1.2 (2) 


— 0.4 


0.18 - 0.22 


-1.20 (8) 


-0.7 (2) 


-0.2 


0.21 - 0.22 


-3.3 (5) 


-0.3 (1) 


-0.1 


0.23 - 0.25 


-5.0 (10) 


-0.2 (1) 



TABLE II: Tricritical points as obtained by the transfer matrix method for several values of H. 
The decoupling point K mi = is included here as the end point of the tricritical line, although it 
does itself not belong to the tricritical three-state Potts universality class. 



H 


K mi 




0.00 


0.0000 (0) 


ln(3)/4 (0) 


0.05 


-0.0107 (12) 


0.269 (1) 


0.10 


-0.0214 (10) 


0.2654 (5) 


0.5 


-0.0937 (5) 


0.2572 (5) 


1.0 


-0.1799 (2) 


0.2500 (2) 


1.5 


-0.2644 (2) 


0.2452 (2) 


2.0 


-0.3481 (2) 


0.2421 (2) 


3.0 


-0.5150 (1) 


0.23845 (8) 


4.0 


-0.6816 (1) 


0.23678 (8) 


5.0 


-0.84823 (5) 


0.23599 (8) 


6.0 


-1.01487 (5) 


0.23560 (8) 
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FIG. 1: The triangular lattice with nearest-neighbor couplings K nn , next-nearest-neighbor cou- 
plings K nnn (examples of which are shown as bold bonds) , and a field H (bold circle) . The lattice 
is divided into three sublattices labeled 1, 2 and 3. 
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FIG. 2: Three lines of phase transitions in the (K nnn , K nn plane. The numerically determined data 
points are shown as circles. The upper line displays the ferromagnetic critical line for K nD > 0. 
For K mi < there are two more lines which represent the boundaries of a critical phase which 
resembles the low-temperature phase of the XY model. The two lines appear to meet at a single 
point, the decoupling point, at K mi = 0. The right hand critical line marks a roughening transition 
to a flat SOS phase, the left hand line a KT-like transition between the disordered and the critical 
phases. The numerical errors in the ferromagnetic region are much smaller than the size of the 
symbols; for the remaining data they are difficult to estimate but believed to be at most of the 
same order as the symbol size. 
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0.12 



FIG. 3: The ferromagnetic critical line, plotted as K nn versus [ln(3)/4 — K nnn ] 7 / 4 . The approximate 
linear behavior confirms that the exponent y a associated with K nn obeys the theoretical prediction 
y a = 7/4. The estimated errors are smaller than the symbol size. 
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FIG. 4: Antiferromagnetic (K nn < 0) critical lines near the decoupling point. The numerical 
results (circles) are plotted as K nn versus [ln(3)/4 — K nnil } 7 / 4 . The approximate linear behavior at 
small |.fC nn | confirms that the exponent associated with the scaling of K mi obeys the theoretical 
prediction y a = 7/4. The estimated errors in the data points are at most of the same order as the 
symbol size. 
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FIG. 5: Finite-size estimates of the magnetic scaling dimension X m versus next-nearest-neighbor 
coupling i^nnn at K nn = — oo. For clarity we include four lines connecting data points for system 
sizes L = 6, 12, 18, 24 respectively. The dashed line indicates the special value X m = 1/18, and 
the black triangle shows the estimated critical value of K nmi for K nn — > — oo. 
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FIG. 6: Finite-size estimates of X m versus if n nn at K nn = —0.6. For clarity we include four lines 
connecting data points for system sizes L = 6, 12, 18, 24 respectively. The dotted and dashed lines 
indicate the special values X m = 1/8 and X m = 1/18 respectively. The two black triangles show 
the estimated critical values of K mm at K mi = —0.6. 
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FIG. 7: Finite-size estimates of the inverse magnetic scaling dimension X~ versus next-nearest- 
neighbor coupling K nnn at K nn = —0.6. The meaning of the lines and symbols are the same as in 
Fig. El The phase transition to flat phase is clearly visible in this figure. 
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nn 

FIG. 8: Dimensionless amplitude ratio Q s versus K nnn at K nn = —0.2. Intersections are found to 
occur near the transition point between the disordered and the algebraic phases. The four lines 
connecting the data points represent, with increasing slope, system sizes L = 24, 36, 48, and 60, 
respectively. The numerical uncertainty margins are much smaller than the size of the data points. 
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FIG. 9: Dimensionless amplitude ratio Q m versus K nnn at K nn = —0.2. Intersections are found 
to occur near the transition point between the algebraic and the flat SOS phases. The four lines 
connecting the data points represent, with increasing slope, system sizes L = 24, 36, 48, and 60, 
respectively. The numerical uncertainty margins are much smaller than the size of the data points. 
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FIG. 10: The complete phase diagram in the three-parameter space (H, K nnn ,e 2Knn ). The solid 
lines denote second-order phase transitions, and the heavy dotted line is the tricritical line sepa- 
rating the three-state Potts critical sheet from the first-order sheet which, shown by heavy dashed 
lines. The three-state Potts critical surface is believed to connect to the e 2Knn = plane at the KT 
line near the origin, and at the K nmi axis until the appearance of the critical phase. The algebraic 
phases for H = and for T = are lightly shaded, and the thin dashed lines are projection lines 
added for clarity. The error margins are at most of the same order as the thickness of the lines. 
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FIG. 11: The tricritical line shown as K Dn versus H. The numerically determined tricritical 
points are shown as circles, and the solid line represents the tricritical lattice-gas limit as K nD = 
-H/6 - 0.01481. 
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